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THE  DEVELOPMENT  OF  AN  ICE-OCEAN  COUPLED  MODEL 
FOR  THE  NORTHERN  HEMISPHERE 


1.0  INTRODUCTION 

The  Polar  Ice  Prediction  System  (PIPS),  based  on  the  Hibler  ice  model  [Hibler  1979;  1980], 
is  an  operational  numerical  model  used  by  the  Fleet  Numerical  Meteorology  and  Oceanography 
Center  (FNMOC)  for  daily  prediction  of  ice-drift  velocity,  thickness,  and  concentration  in  the 
Arctic  Ocean.  PIPS  uses  daily  information  on  ocean  currents,  oceanic  heat  flux,  incoming  solar 
radiation,  sensible  heat  flux,  total  heat  flux,  surface  air  temperatures,  geostrophic  wind  stresses,  and 
new  data  on  ice  concentrations  to  make  these  estimates.  The  oceanic  variables  that  it  needs,  ocean 
currents  and  heat  fluxes,  are  monthly  mean  values  derived  from  the  Cox  ocean  model  [Hibler  and 
Bryan  1987]. 

The  PIPS  domain  mainly  covers  the  central  Arctic,  the  Barents  Sea,  the  East  Greenland  Sea, 
and  the  Norwegian  Sea  (Fig.  1.1).  The  PIPS  grid,  a  subsection  of  the  FNMOC  polar  stereographic 
projection  grid  for  the  Northern  Hemisphere,  has  a  dimension  of  47  x  25  grid  squares,  each 
approximately  127  km  on  a  side.  Since  the  difference  between  the  spherical  grid  and  the  stereo¬ 
graphic  projected  grid  is  minimal  in  the  Arctic,  less  than  0.8%  throughout  the  model  domain,  the 
PIPS  model  Cartesian  grid  can  be  used  to  estimate  sea  ice  conditions  without  losing  much  accuracy. 

South  of  the  Arctic,  however,  in  areas  such  as  the  Labrador  Sea  and  the  Sea  of  Okhotsk,  the 
differences  between  the  spherical  grid  and  stereographic  projected  grids  are  significant  enough  that 
the  PIPS  model  cannot  be  simply  extended  to  accurately  compute  seasonal  conditions  there.  The 
PIPS  model  must  be  transformed  into  spherical  coordinates,  and  the  most  appropriate  transformation 
is  discussed  in  this  report. 

Owens  et  al.  [1990]  transformed  the  Hibler  ice  model  to  the  spherical  coordinates  of  the  Earth 
latitudes  and  longitudes  for  Antarctic  sea  ice.  Since  the  land  area  of  the  Antarctic  continent  covers 
most  areas  from  80°  S  to  the  South  Pole,  the  Owens  et  al.  ice  model  had  no  numerical  singularity 
at  the  Pole,  and  very  little  numerical  instability  in  high  latitudes.  Thus,  there  was  no  need  to  solve 
the  flow  equations  in  the  actual  polar  areas  in  the  Owens  et  al.  model. 

However,  this  set  of  spherical  coordinates  would  not  work  for  the  Northern  Hemisphere  because 
of  the  circulation  calculations  required  for  the  North  Pole.  Flato  and  Hibler  [1993]  and  Holland 
[1994]  have  tried  to  use  the  standard  Earth  latitudes  and  longitudes  in  models  of  the  Arctic.  They 
introduced  an  “island”  at  the  North  Pole  to  avoid  making  calculations  there,  but  the  results  show 
sea  ice  obviously  drifting  around  the  island  rather  than  through  it.  Holland  then  transformed  the 
North  Pole  to  Siberia  for  calculations,  and  sea  ice  drift  across  the  Pole  occurred.  Since  they  used 
the  same  atmospheric  forcing  for  both  sets  of  calculations,  the  difference  must  have  been  caused 
by  the  coordinate  systems  that  were  used. 

In  addition  to  the  problem  of  sea  ice  not  drifting  across  the  Pole,  sea  ice  was  thinner  in  the 
island  models  than  in  the  transformed  pole  model  calculations.  The  sea  ice  was  up  to  0.5  m  thicker 
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in  the  transformed  pole  model  calculations  than  in  the  island  model  calculations.  The  distribution 
of  ice  thickness  was  similar  in  the  transformed  Pole  model  calculations  to  other  work  [Cheng  and 
Preller  1992;  1994]:  ice  was  thick  north  of  Greenland  and  the  Canadian  Archipelago,  less  thick 
across  the  North  Pole,  and  thin  near  the  Russian  coast.  Consequently,  for  the  PIPS  model  to 
describe  the  evolution  of  conditions  throughout  the  Northern  Hemisphere  it,  too,  needs  the  transformed 
pole  of  the  spherical  coordinate  system  rotated  away  from  the  geographic  North  Pole. 

This  report  describes  the  transformation  of  the  PIPS  model  from  rectangular  coordinates  to  the 
new  spherical  coordinates,  and  the  PIPS  model  modifications  required  to  accommodate  that  trans¬ 
formation.  The  transformed  model  grid  matches  the  original  PIPS  grid  closely  where  they  overlap, 
and  it  can  estimate  the  evolution  of  ice  conditions  south  of  the  original  PIPS  grid  without  the 
numerical  difficulties.  In  the  following  discussions,  we  first  show  the  transformation  that  was  made, 
and  then  test  cases  that  compared  the  results  of  the  new  model  against  the  rectangular  PIPS.  It  is 
noted  how  the  transformed  model  maintains  the  same  physics  as  the  rectangular  PIPS,  and  the 
programming  changes  necessary  to  maintain  the  physics  are  described.  Test  results  are  divided  into 
the  model  domain’s  seven  subregions:  the  Sea  of  Okhotsk,  the  Bering  Sea,  the  central  Arctic,  the 
Barents  Sea,  Hudson  Bay,  the  Labrador  Sea/Baffin  Bay,  and  the  Norwegian/East  Greenland  Seas. 
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The  transformed  PIPS  model  was  extended  southward  (to  20°  N-30°  N,  depending  on  location) 
to  cover  most  sea  ice  formation  and  drift  in  the  Northern  Hemisphere  with  fine  grid  resolution.  The 
coupling  of  the  ice  model  to  the  Cox  ocean  circulation  model  for  calculation  purposes  is  described, 
as  well  as  how  the  atmospheric  forcing  fields,  bathymetry  data,  ocean  temperature,  and  salinity 
fields  were  interpolated  to  support  studies  with  the  transformed  and  extended  PIPS  model  at  lower 
latitudes.  The  autonomy  of  the  PIPS  model  and  the  ocean  circulation  model  (the  Cox  model  [1984] 
for  this  work)  was  maintained  during  this  development  as  much  as  possible  so  that  the  Cox  ocean 
model  could  be  replaced  in  future  developments  by  upgraded  code  (be  it  an  upgraded  Cox  model, 
or  the  Mellor  ocean  model  [1992]). 

At  present,  there  are  three  operational  PIPS  model  versions:  a  coarse  spatial  grid  version  (127  km), 
which  covers  the  entire  Arctic  area  of  interest,  and  two  finer  grid  resolution  versions  (20-km 
grid  squares):  one  for  the  Barents  Sea  (the  Regional  Polar  Ice  Prediction  System-Barents  Sea 
[RPIPS-B])  and  one  for  the  East  Greenland  Sea  (the  Regional  Polar  Ice  Prediction  System-Greenland 
Sea  [EPIPS-G]).  The  finer  grid  models  are  not  coupled  with  ocean  models;  they  are  driven  by 
constant  ocean  currents  and  oceanic  heat  fluxes,  and  will  be  replaced  in  the  future  by  a  finer  grid 
ice-ocean  coupled  model. 

The  ice-ocean  coupled  model  has  been  applied  to  studies  other  than  operational  uses,  such  as 
modeling  the  changes  in  the  Arctic  Ocean  due  to  variable  river  outflows  [Allard  et  al.  1995].  In 
addition,  the  “ocean-only”  part  of  the  model  has  been  used  to  study  the  dispersion  of  radioactive 
waste  material  in  the  Kara  and  Barents  Seas  in  the  Arctic  and  in  the  North  Atlantic  [Preller  et  al. 
1993;  Preller  and  Cheng  1995].  It  is  anticipated  that  further  applications  of  the  model  will  evolve 
with  time. 


2.0  NEW  SPHERICAL  COORDINATES 

Figures  2. 1  and  2.2  show  a  coordinate  transformation  from  the  Earth  latitudes  and  longitudes 
to  a  new  spherical  coordinate  system  that  removes  computational  instabilities  from  the  PIPS  model. 
This  particular  coordinate  transformation  has  the  following  advantages.  There  is  neither  numerical 
instability  in  high  latitudes  nor  numerical  singularity  at  the  North  Pole  of  this  coordinate  system 
because  the  new  north  pole  has  been  transformed  to  near  the  equator.  In  addition,  putting  the  new 
north  pole  in  a  land  mass  means  computations  near  or  across  a  pole  will  never  be  an  issue.  For  the 
PIPS  domain,  the  spherical  grid  to  which  the  computational  algorithms  were  transformed  and  the 
Cartesian  grid  of  the  original  formulation  are  almost  identical.  Therefore,  when  the  extra  terms  are 
added  to  account  for  the  coordinate  transformation,  the  physics  remain  the  same,  and  one  grid  can 
be  verified  against  the  other.  The  same  spherical  geometry  was  introduced  into  the  Cox  ocean 
model  [1984],  then  both  the  PIPS  ice  model  and  the  Cox  ocean  model  were  coupled  in  spherical 
coordinates. 

The  new  set  of  spherical  coordinates  for  the  PIPS  model  was  constructed  in  the  following 
manner.  First,  longitudes  were  rotated  so  that  the  prime  meridian  was  located  at  190°  E  (Fig.  2.1). 
Then,  the  North  Pole  was  rotated  90°  down  the  100°  E  meridian  until  it  resided  on  the  true  equator 
(Fig.  2.2).  Therefore,  the  transformed  coordinates  coincide  with  the  170°  W-10°  E  great  circle, 
which  passes  through  the  North  Pole  as  the  new  equator.  The  new  equator  passes  right  through  the 
old  PIPS  grid,  where  each  grid  square  in  the  new  spherical  coordinates  is  about  1.143°  on  a  side. 
The  Cartesian  PIPS  grid  and  the  spherical  coordinates  PIPS  grid  are  almost  identical.  The  X-  and 
Y-axes  are  stereographic  projections  of  the  170°  W— 10°  E  and  the  80°  W— 100°  E  great  circles, 
respectively. 
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Fig.  2.1  — Coordinate  transformation  from  the  Earth-oriented  longitudes  and  latitudes, 
step  one:  rotate  Greenwich  meridian  190°  eastward 


To  convert  the  model  from  Cartesian  coordinates  to  the  new  spherical  coordinates,  replace 
and  63:  by  6(t>  and  60  as  follows: 

8x  =  r  cos06(t)  > 

By  =  r  50  , 

where  Bx  and  6y  represent  spatial  intervals  in  Cartesian  coordinates,  r  is  the  Earth’s  radius,  0  is  the 
latitude  in  the  spherical  coordinates,  (j)  is  the  longitude,  and  6<1)  and  60  are  increments  in  the  spherical 
coordinates.  The  physics  of  the  ice  and  ocean  models  were  recast  into  these  variables. 

The  edges  of  the  PIPS  grid  are  about  1 1  grid  squares  from  the  new  equator  of  the  transformed 
coordinates,  or  about  12.57°.  The  difference  between  the  Cartesian  grid  (a  projected  grid)  and  the 
spherical  grid  is  less  than  0.8%  at  the  boundaries  (i.e.,  the  difference  between  the  12.57°  arc 
distance  and  the  sin  12.57°).  Section  5.0  of  this  report  provides  a  more  thorough  comparison 
between  the  Cartesian  model  and  the  model  with  the  transformed  spherical  coordinates.  The  results 
of  ice-drift  velocity,  concentration,  and  thickness  cannot  be  distinguished  in  vector  or  contour  plots. 

3.0  ICE  MOMENTUM  BALANCE  EQUATION 

As  in  Hibler  [1979],  Preller  and  Posey  [1985],  Hibler  and  Bryan  [1987],  and  Cheng  and  Preller 
[1992;  1994],  the  momentum  balance  equation  for  sea  ice  is  as  follows: 
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Fig.  2.2  —  Coordinate  transformation  from  the  Earth-oriented  longitudes  and 
latitudes,  step  two:  rotate  North  Pole  90°  southward  along  the  100°  E  meridian 


3u 

+  mu«Vu  +  w/kxu  =  , 

at 


(1) 


where  m  represents  ice  mass  per  unit  area, 

u  is  ice-drift  velocity  expressed  in  spherical  coordinates  as 

$  and  6  are  the  unit  vectors  in  the  angular  coordinates  of  the  spherical  coordinate  system  into 
which  the  equation  of  motion  must  be  transformed, 

/  is  the  Coriolis  parameter, 

k  is  the  unit  vector  that  points  radially, 

Xa  is  the  wind  stress  at  the  ice  surface, 

Xw  is  the  ocean  current  drag  stress  on  the  bottom  of  the  sea  ice, 
g  is  the  gravitational  acceleration. 
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H  is  the  dynamic  height  of  ocean  currents,  and 
F  is  a  force  produced  by  ice  internal  stress. 
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This  equation  must  be  transformed  into  spherical  coordinates  to  ensure  that  the  physics  of  the 
ice-ocean  model  is  correct. 


du 

As  long  as  u  is  expressed  in  spherical  coordinates,  the  temporal  derivativem term  needs  no 
modification  for  inclusion  in  the  PIPS  model. 

The  gradient  term  V  is  modified  in  the  following  way.  Since  the  sea  ice  drifts  horizontally  at 
the  sea  surface,  the  vertical  derivative  is  zero  in  the  gradient  term.  The  other  components  are 


wntten  as 


f  1  ^  gl  — 

^  r  COS0  3(t)  ’  rdB 


.  The  total  advection  term  is  equivalent  to  mV  up 


ux ( V  XU ) 


in  spherical  coordinates  where  the  x  is  the  cross-product.  Expanding  this  term  gives,  for  the  advection 


f 


term,  m 


1  3u  1  du  ^UQU^tand  M^uxtane 

- -  ^  +  Ua  -  —  -  6 - - +  9 

r  cos6  ^  r  dQ  r  r 


\ 


The  Coriolis  term  m/kxu  is  written 


) 


as  mf  (kxgw^-i-  kx$Me).  Note  that  the  Coriolis  parameter  should  be  evaluated  for  each  grid  cell 
because  of  the  variations  in  its  value  with  latitude. 

The  terms  for  wind  and  ocean  current  stresses  are  expressed  as  follows: 
t„=Pa,vC.|U,|(u  g  COS  a  ^  +  kxU^  sina  and 


^w~  Pw^v 


where  represents  the  geostrophic  winds, 

Uh,  is  the  ocean  current  of  the  mixed  layer. 

Pair  and  Piv  are  the  densities  of  air  and  seawater, 

Ca  and  Cw  are  drag  coefficients  of  wind  and  ocean  currents,  and 

tta  and  aw  are  the  turning  angles  between  air  and  ice  and  between  ice  and  water,  respectively. 

The  law  of  the  wall  was  applied  to  compute  the  variable  coefficient  Cw  as  a  function  of  sea  ice 
thickness.  Further  discussions  and  comparisons  are  presented  in  Secs.  12.0  and  14.0.  In  those 
discussions,  the  symbol  Cdjce-water  is  used  for  the  drag  between  sea  ice  and  water. 

The  force  F  due  to  the  ice  internal  stress  aij  can  be  written  as  follows: 


1  a 

rcosG  90 


(COS0a0y)  + 


1 

rcosG  3(t) 


or 
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F  =  V.a,.,-  =  ( 


-  tan0  1 

— — + 
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da. 


r  90  rcos0  9(t) 


-)$  +  (■ 


tan0 


<^69  + 


1  9o 


00 


1 


9<y. 


r  90  rcos0  9(t) 


-)§,  (2) 


where  i  and  j  represent  the  indices  of  the  variables  0  and  (j).  Note  that  the  divergence  operator 

introduces  the  extra  terms  and  ^^^^04,  when  it  operates  due  to  the  coordinate 

r  r  ^ 

transformation.  The  ice  internal  stress  (])//  and  the  corresponding  strain  £y  (or  a  strain  rate  to  be 
precise)  are  shown  below; 

Oij  =  2riey  +  5y  ((^  -  n)  eu-  -  ■P/2), 


where  i;  =  P/2A  and  ri  =  ^le'^  (computed  from  and  u%  in  the  previous  time  step). 


i 


1/2 


+  2  £aa£ 


00 


P  =  2.75« 

A  and  h  are  ice  concentration  and  thickness,  respectively, 

8y  is  the  Kronecker  delta  function:  1  as  i  =j  and  0  otherwise, 


e  =  2. 


^kk  ~  ^00  ^00  ’ 


1  W0tan0 


rcosB  00 


^  00  “  ^  60  ““  2^ 


1  9m  e  M  AtanO 

+  — - +  — ,  and 


’I 


'COS0  9(1) 


r 


2  _  1  0 
^00“7'^- 

do 

Therefore,  the  force  — ^  due  to  the  ice  internal  stress  can  be  expressed  as  follows  (6x,-  and  bxj 
dxi 

represent  either  rcos05(l)  or  r56); 
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1  un  tan0 

■(2il(— -77-rr  +  — )) 


rcos09(t)  rcos0  0(t) 


a  1  1  M<|)tane  1  du^ 


3  1  dufy  Metan0  2  duQ 

*  7^  « ^  - -I '  <  7^59  ^  " -r- ^  7  M  > 

a  .  1  a  1 

rcosSac])  rcos9  a(t)  ^  ra0  ^  rcos0  a(j) 


+  ^  (11  -  ^ )+ c- n )- ^) 


r30  r  30  rcos03(j) 


rcos0a(t) 


r  rcos0a(j) 


r  30 

(2ri«e) 


tan0  3(tiha)  sec^0  tan0  d  ^ 

+ - ^  + - ^  + - E3l((^-il)“e)- 

r  rod  r  r  r  r cos 03(1) 


Following  Hibler  [1979],  apply  the  finite  difference  method  to  a  discretized  form  of 
equation  above  becomes: 


3;C; 


^®i<j)  1 

dxj  2 


*^(j)m-In(ll/nn  +  l  '^dm.n'^^mn  +  X 
1  j  ~  (|)/nn  (Hm/?  +  I"*"!!  mn  1  «  +  1  ■*■11 /n+ I  « 

rcos0A(^)  \  +  i  +  ^ 

mn  ^/n  +  ln  +  l'^Cm+lw) 

^  +  1  m+1/2  +  1  '^^m+ln'^^m+in  +  l  +  I  «) 


1 

+  “ 


)^m+ 1  w  +  1  0m+ 1 /2  +  1  ^  0m  + 1 /2  ^Bmn  +  l  ^  %mr^ 

m  /z  +  1  ^  dm  /I  +  1*^  ^  0An  n~  ^  dm  -  \  n  +  l~  ^  dm~\ 
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l/  1  vz  ^  ^  ^  ^ 

^1  7A0  )  \  ^  ^  n^m  Az“*"^m+lAz‘^^m  Az  +  l'^^m+lAz  +  l) 
“*"^cl)mAz  +  l(nmAz  +  i  az+i) 


and  the 


The  Development  of  an  Ice-Ocean  Coupled  Model 


9 


1 

+  — 


1 


4  .2 


r^cos  0A(t)A0 


1 


(Cm 

+  lrt 

+  1 

-Tl 

+  (C 

m+  1 

/z 

-Tli 

+  (C 

m  /?  + 1 

+  (C 

m/i 

n 

m/z) 

tan0 

1 

rcos  0A(}) 


1‘^Qm  +  ln(^m  n  +  l'^^m  n^~  +  -  l). 


H - 4  I  ^  9'””  +  1  (^m+  1  n  ^  t 


sec^0  1  ,  ^  \ 

+  — ~ —  ^  W  ^mn  Ol/nn  "*■  ^  mn  +  1  ^  m+  1  n  ■*"  ^  m+  1  n  +  I-' 


tan0  1  1  j  ^  6m+ I  n  (^mn  +  I  ^mn  ^mn+1  ^mn) 

r  rCOS0A(j)  4  1  „(Cm«  +  Cm  n  -  1  “  "nmn  ~  'H/nn-  l)  | 


+ 


Similarly,  the  force 


9j:, 


due  to  the  ice  internal  stress  can  be  written  as  follows: 
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The  discretized  form  of 


above  becomes: 
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4.0  THE  CONTINUITY  EQUATION  FOR  ICE  CONCENTRATION  AND  THICKNESS 

The  continuity  equation  for  ice  concentration  is  written  as  follows  [Hibler  1979]: 

—  +  V«  (uA]  =  +  (diffusion),  (3) 


where  A  represents  the  ice  concentration, 

V»(uA)  is  the  advection  term  for  ice  concentration, 

Sa  is  the  ice  growth  rate  (i.e.,  a  thermodynamic  term  due  to  atmospheric  and  oceanic  forcing), 
and  the  diffusion  term,  usually  small,  is  mainly  used  to  stabilize  the  numerical  algorithm. 

dA 

Again,  as  with  Eq.  (1),  the  temporal  term  would  not  be  affected  by  transforming  the 
coordinates  from  Cartesian  to  spherical.  Its  form  remains  the  same  in  the  spherical  version  of 
the  continuity  equation. 

The  advection  term  V»(uA)  in  spherical  coordinates  is  transformed  as  follows: 


„  -tane  ,1  ^(«e^)  1 

V.(u;l)  =  -— +  3<|>  ' 


The  extra  term  — ^^^M9Ais  a  result  of  the  coordinate  transformation. 


The  diffusion  term  is  a  sum  of  two  terms  V^A  and  V^V^A,  which  transform  as  follows: 


-tane  dA  1  3^A  1  9^A 

^2  90  ^2  002  r^cos^Q  3<t)^ 


30^  r^cos^0  3(|)^ 

The  thermodynamic  term  Sa  is  expressed  as  follows: 

Sa  =  ^-^  (1- A)//(  /(O))  +  ^ShHifih/A)), 

rtQ 

where  /(O)  is  the  growth  rate  of  ice  thickness  in  open  water, 

H(x)  is  the  Heaviside  step  function,  and  equals  1  for  x  >  0  and  0  otherwise, 
h  represents  the  ice  thickness,  and  ho  is  set  to  be  0.5  m, 

Sh  is  the  net  ice  growth  or  melt,  and  equals  f{h/A)A  +  /(0)(1  -  A),  and 
f{h/A)  is  the  growth  rate  of  ice  thickness  h/A. 


12 


Cheng  and  Preller 


The  function  /(O)  becomes  the  growth  rate  due  to  the  atmospheric  cooling  after  all  the 
mixed-layer  heat  above  the  freezing  temperature  is  used  [Cheng  and  Preller  1992;  1994],  The  mixed- 
layer  temperature  used  here  comes  from  the  Cox  ocean  model  that  is  the  basis  of  the  ice-ocean 
coupled  model. 

Note  that  the  Hibler  ice  model  [Hibler  1979;  1980]  calculates  the  atmospheric  heating  of  the 
mixed  layer  as  a  negative  ice  growth  rate.  In  PIPS2.0,  that  heating  comes  to  the  mixed  layer  from 
the  ocean  model.  Further  discussions  about  the  negative  ice  growth  rate,  the  mixed-layer  temperature, 
and  the  oceanic  heat  flux  can  be  found  in  Sec.  12.0  of  this  report. 

All  changes  discussed  above  for  the  continuity  equation  for  ice  concentration  apply  as  well  to 
the  continuity  equation  for  ice  thickness.  The  derivation  of  the  thickness  equation,  which  has  been 
omitted  here,  is  completely  analogous,  except  that  the  Sa  term  is  replaced  by  S/j  for  the  thickness 
equation. 


5.0  MODEL  COMPARISON  IN  THE  PIPS  DOMAIN 

The  spherical  coordinate  version  of  PIPS  was  compared  here  with  the  Cartesian  version  of  PIPS 
to  ensure  that  no  discernible  inaccuracies  were  introduced  during  the  coordinate  transformation. 
During  the  test,  1986  NOGAPS  winds  were  used  for  atmospheric  forcing  and  the  monthly  oceanic 
forcing  was  provided  by  Hibler  and  Bryan  [1987].  Both  the  Cartesian  and  the  spherical  coordinate 
ice  models  were  spun  up  for  3  yr  with  the  following  initial  conditions:  3-m  ice  thickness,  100%  ice 
concentration,  and  zero  ice  velocity.  Figures  5.1a,  5.2a,  and  5.3a  show  ice  velocity,  thickness,  and 
concentration,  respectively,  for  17  Jan  1986,  computed  from  the  spherical  coordinate  PIPS  model. 
Figures  5.1b,  5.2b,  and  5.3b  show  the  same  quantities  for  the  Cartesian  model.  There  are  no  visible 
differences  between  Figs.  5.1a  and  5.1b,  5.2a  and  5.2b,  and  5.3a  and  5.3b,  although  numerical 
differences  of  a  few  percent  exist.  The  model  comparison  demonstrates  that  (1)  both  models  give 
almost  identical  results,  and  (2)  the  transformed  PIPS  model  can  be  applied  to  cover  most  sea  ice 
regions  in  the  Northern  Hemisphere  without  numerical  difficulties  [Cheng  and  Preller  1992; 
Holland  et  al.  1994]. 


6.0  PIPS  DOMAIN  IN  THE  NORTHERN  HEMISPHERE 

Recasting  the  ice-ocean  model  in  spherical  coordinates  provides  a  different  polar  grid  on  which 
to  calculate  ice  and  ocean  conditions  from  the  one  that  was  used  for  the  Cartesian  model  calculations 
(see  Fig.  1.1).  Figure  6.1  shows  the  PIPS2.0  model  domain  on  a  stereographic  projection  of  the 
Northern  Hemisphere  (the  hatched  area).  The  domain  includes  the  Gulf  of  St.  Lawrence,  the  Labrador 
Sea,  Baffin  Bay,  Hudson  Bay,  the  Baltic  Sea,  the  Yellow  Sea,  the  Sea  of  Okhotsk,  the  Bering  Sea, 
and  many  other  seas  of  the  Arctic  and  sub-Arctic  (Fig.  6.2).  This  domain  allows  no  sea  ice  inflow 
or  outflow  as  boundary  conditions  because  no  sea  ice  exists  as  far  south  as  20°-30°  N.  (Note  that 
since  sea  ice  is  the  main  interest,  the  following  seas  and  lakes  are  ignored:  the  five  Great  Lakes 
of  North  America,  the  Mediterranean  Sea,  and  the  Caspian  Sea.) 

Figure  6.1  was  drawn  for  eveiy  fourth  grid  line  and  has  360  x  360  grid  cells,  each  with  an  arc 
length  of  approximately  0.28575°  for  both  sides  of  the  grid  square.  Note  that  at  the  North  Pole,  the 
equal  angular  measures  for  the  grid  square  sides  translate  into  equal  linear  dimensions  of  approximately 
32  km.  However,  the  grid  cell  at  the  top  boundary  of  Fig.  6.1  shows  the  most  linear  inequality  for 
the  grid  dimensions  (i.e.,  approximately  17  x  32  km).  In  the  new  model  domain,  the  equator  that 
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Modeled  ice  concentration  from  PIPS  model  coordinates:  (a)  spherical 
and  (b)  Cartesian.  Concentration  ranges  are  from  0  to  1. 
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Fig.  6,1  —  PIPS2.0  domain  with  the  variable  resolution  grid 
overlaid.  The  hatched  lines  are  drawn  at  every  fourth 
grid  point  including  land  points. 


Fig.  6.2  —  The  PIPS2.0  model  domain,  is  equivalent  to  the  hatched  area  in  Fig.  6.1, 
is  plotted  with  the  same  arc  distance,  0.28575°,  for  both  vertical  and  horizontal  axes. 
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is  defined  by  the  polar  rotation  lies  on  the  170°  W— 10°  E  great  circle  that  passes  through  the 
North  Pole. 


7.0  ATMOSPHERIC  FORCING 

The  atmospheric  forcing  includes  incoming  solar  radiation,  sensible  heat  flux,  total  heat  flux, 
surface  air  temperature,  atmospheric  pressure,  and  vapor  pressure.  Oceanic  and  atmospheric  heat 
fluxes  contribute  to  a  heat  balance  equation  that  allows  the  calculation  of  ice  growth/decay.  Ocean 
currents  in  the  top  mixed  layer  and  the  geostrophic  wind  stress  are  used  to  calculate  ice  motion. 

The  1986  daily  NOGAPS  atmospheric  forcing  was  used  to  compute  monthly  averages  that  were 
interpolated  into  the  P1PS2.0  domain  using  two-dimensional  cubic  splines  to  provide  the  necessary 
forcing  fields.  As  examples  of  the  kinds  of  data  produced  by  this  procedure.  Figs.  7.1  and  7.2  show 
the  monthly  NOGAPS  geostrophic  winds  and  surface  air  temperature  averages,  respectively.  The 
winds  were  computed  from  the  pressure  fields  using  the  geostrophic  approximation.  The  wind 
velocities  were  plotted  at  every  fourth  point  of  the  P1PS2.0  grid. 

Because  of  such  weather  features  as  cyclones,  wind  magnitudes  and  directions  vary  throughout 
the  year  in  the  Arctic.  Sea  ice  motion  changes  follow  accordingly.  Westerly  winds  prevail  for  most 
of  the  year  in  the  North  Pacific  and  North  Atlantic  Oceans,  but  the  winds  become  somewhat  more 
random  during  the  summer. 

Long-wavelength  radiation  was  calculated  based  on  the  total  heat  flux,  the  sensible  heat  flux, 
and  the  incoming  solar  radiation  from  the  NOGAPS  data  fields.  Boltzmann  radiation  is  included  in 
the  computed  long-wavelength  radiation.  The  incoming  solar  radiation  is  defined  as  positive  in  all 
calculations  (i.e.,  heat  flow  into  the  sea  ice  and  water  is  positive).  Small  negative  values  that  occur 
during  the  balance  calculations  with  the  other  heat  sources  are  probably  caused  by  errors  in  the 
interpolation  from  the  NOGAPS  model  grid  to  the  ice/ocean  model  grid.  Therefore,  these  small 
errors  are  set  to  zero  in  the  calculations. 

General  comments  describing  the  Arctic  input  data  will  be  useful  in  understanding  the  results 
of  running  the  sea  ice  model  that  follows.  Incoming  solar  radiation  decreases  with  latitude.  In  the 
Arctic,  vapor  pressures  are  small  from  late  fall  to  early  spring,  but  reach  up  to  10  mbar  during 
the  rest  of  the  year.  In  general,  vapor  pressures  decrease  with  latitude.  Sensible  heat  fluxes  seem 
to  be  randomly  distributed  in  the  Arctic  throughout  the  year,  but  they  generally  increase  south  of 
the  Arctic. 

The  NOGAPS  forcing  fields  in  the  63  x  63  FNMOC  grid  of  the  Northern  Hemisphere  have  a 
horizontal  resolution  for  each  grid  cell  of  approximately  380  km.  That  spatial  resolution  is  too 
coarse  to  resolve  some  regions  (such  as  the  Baltic  Sea,  the  Yellow  Sea,  the  Sea  of  Okhotsk,  the 
Gulf  of  St.  Lawrence,  the  Labrador  Sea,  and  Baffin  Bay).  The  NOGAPS  data  must  be  interpolated, 
which  can  introduce  error  into  the  fields  for  grid  squares  close  to  land.  Estimates  of  ice  drift 
velocity,  thickness,  and  concentration  produced  by  the  research  version  of  PIPS2.0,  which  includes 
the  interpolation,  could  be  biased.  In  operation,  when  implemented  at  FNMOC,  the  coupled  model 
(PIPS2.0)  will  use  the  NOGAPS  fine  spatial  resolution  data  fields  (i.e.,  1°  x  1°  in  longitudes  and 
latitudes)  for  calculations.  The  error  in  interpolated  forcing  fields  should  be  reduced  considerably. 

Another  source  of  error  in  the  PIPS2.0  model  has  to  do  with  how  the  wind  stress  is  implemented. 
Instead  of  surface  winds  and  their  stresses,  monthly  geostrophic  winds  and  a  corresponding  drag 


Fig.  7.1  — Monthly  NOGAPS  geostrophic  winds  (m/s),  Jan-Dec  1986 


Fig.  7.1  — cont. 
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Fig.  7.1 — com. 
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Fig.  7.2  —  Monthly  NOGAPS  surface  air  temperature,  Jan-Dec  1986  (°C) 
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Fig.  7.2  —  cont. 
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Fig.  7.2  —  cont. 
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coefficient  0.001  were  used  to  drive  the  coupled  model  as  an  expedient  way  of  speeding  up  the 
model  execution  and  as  a  useful  and  accurate  approximation  for  some  ocean  areas.  In  the  Arctic, 
the  North  Pacific,  and  the  North  Atlantic  Oceans,  this  geostrophic  approximation  works  reasonably 
well  because  of  the  expanses  of  open  water.  However,  the  accuracy  of  the  results  is  suspect  near 
coasts,  such  as  the  southeast  Bering  Sea  between  Alaska  and  the  Aleutian  Islands,  because  the  wind 
estimates  are  affected  by  the  water/land  interfaces.  When  implementation  of  the  PIPS2.0  model  is 
complete  at  FNMOC,  the  geostrophic  wind  and  pressure  fields  will  be  replaced  by  NOGAPS 
surface  wind  stress,  which  will  address  this  problem  for  operations. 


8.0  THE  COX  OCEAN  MODEL 

In  his  original  work  on  Arctic  ice  modeling  [1987],  Hibler  used  the  Cox  ocean  model  [1984] 
to  provide  the  oceanic  forcing  for  the  water  components  of  the  momentum  and  heat  flow  calcula¬ 
tions.  The  PIPS  ice  growth  model  retained  this  dependence  on  the  Cox  ocean  model  although 
mainly  for  the  sake  of  convenience,  as  the  PIPS  investigators  were  occupied  with  including  other 
features  in  the  model.  The  Cox  ocean  model  provides  oceanic  forcing  for  the  ice  models;  namely, 
the  ocean  currents,  the  mixed-layer  temperature  and  salinity,  and  the  oceanic  heat  fluxes.  This 
section  describes  the  version  of  the  Cox  ocean  model  that  was  used  for  PIPS2.0,  as  well  as  the 
procedures  that  were  used  to  execute  it. 

The  Cox  ocean  model  is  a  three-dimensional,  primitive  equation  numerical  model  for 
large-scale  baroclinic  ocean  circulation.  The  model  was  spun  up  from  a  motionless  initial  ocean 
state  for  5  yr  of  model  time  to  develop  steady-state  oceanic  characteristics.  These  characteristics 
were  used  as  initial  values  for  the  ice-ocean  coupled  model.  Sarmiento  and  Bryan’s  robust  method 
[1982]  was  used  during  the  spin-up  to  constrain  temperature  and  salinity  in  the  top  mixed  layer  with 
30  day,  and  then  to  constrain  temperature  and  salinity  in  the  rest  of  the  levels  with  3  yr. 

The  parameters  and  data  for  the  spin-up  initialization  are  as  follows.  The  coefficients  of  horizontal 
eddy  diffusion,  vertical  eddy  diffusion,  horizontal  eddy  viscosity,  and  horizontal  eddy  viscosity  are 
10^cm^s“^  1  cm^s~^  10^  cm^s“^  and  1  cm^s”^  respectively.  The  timestep  is  1  hr  for  temperature 
and  salinity,  and  0.1  hr  for  ocean  currents  and  stream  functions  (i.e.,  the  distorted  physics  method 
of  Bryan  [1984]).  The  1986  annual  NOGAPS  geostrophic  winds  and  the  Levitus  annual  climatological 
data  [1982]  were  used  in  this  5-yr  spin-up  process. 

Assumptions  made  in  the  implementation  of  the  Cox  model  manifest  themselves  with  signatures 
in  the  data  that  the  PIPS  model  produces  during  runs.  Those  assumptions,  as  well  as  the  signatures  that 
they  produce,  are  discussed  in  Secs.  13.0  and  14.0. 


9.0  OCEAN-BOTTOM  TOPOGRAPHY 

To  support  studies  using  the  Cox  ocean  model,  the  NOAA-NGDC  Earth  Topography  System, 
Version  5  (ETOP05)  bathymetric  data  were  interpolated  into  the  PIPS2.0  domain  and  grid  using 
the  cubic  splines.  Figure  9.1  shows  the  ocean-bottom  topography  of  this  model  domain.  The  topog¬ 
raphy  was  smoothed  to  avoid  numerical  difficulties  in  the  Cox  ocean  model  by  integrating  the  cubic 
spline  function  and  then  averaging  it  over  every  3x3  grid  squares  to  represent  the  center  square’s 
topography. 
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The  topography  was  divided  into  15  levels,  according  to  the  requirements  of  the  Cox  ocean 
model  [1984],  Several  channels  and  straits  were  too  narrow  to  appear  as  waterways  after  this 
averaging  technique  was  used.  Examples  are  the  English  Channel  and  the  channel  between  the 
fresh-water  Baltic  Sea  and  the  salt-water  North  Sea.  In  addition,  some  small  islands  were  averaged 
with  adjacent  deep  trenches  and  became  sea  gridpoints.  Where  it  is  important  to  have  land  topography, 
island  points  such  as  some  of  the  Kuril  Islands  and  the  Aleutian  Islands,  have  been  manually  edited 
back.  Table  9.1  lists  the  level  depths  and  the  thickness  of  each  ocean  layer  for  bottom  topography 
used  by  the  Cox  ocean  model.  Any  bottom  topography  from  the  data  set  deeper  than  the  Level  15 
definition  (5700  m)  was  set  to  that  depth  for  this  study.  Similar  adjustments  were  made  on  the 
Levitus  climatological  data  of  temperature  and  salinity  (see  Sec.  10.0)  when  applied  to  studies 
involving  the  PIPS2.0  domain.  Model  ocean  currents  at  depth  could  be  biased  by  such  adjustments 
near  and  in  deep  trenches. 
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Table  9.1  —  Level  Depths  and  Thicknesses 
of  Bottom  Topography  for  the  Cox  Ocean 
Model 


Level 

Thickness 

(m) 

Center  Depth 
(m) 

Depth 

(m) 

1 

30.0 

15.0 

30.0 

2 

46.3 

53.1 

76.3 

3 

67.4 

110.0 

143.7 

4 

94.5 

191.0 

238.2 

5 

128.5 

302.4 

366.7 

6 

170.6 

451.9 

537.2 

7 

221.9 

648.2 

759.1 

8 

283.5 

900.0 

1042.7 

9 

356.2 

1220.8 

1398.9 

10 

440.6 

1619.1 

1839.4 

11 

536.8 

2107.8 

2376.2 

12 

644.5 

2698.5 

3020.8 

13 

763.0 

3402.3 

3783.8 

14 

890.7 

4229.1 

4674.5 

15 

1025.5 

5187.2 

5700.0 

Level  1,  a  30-m-thick  layer  that  interacts  with 
sea  ice  and  atmosphere,  was  chosen  to  represent  an 
average  depth  for  the  mixed  layer  in  the  Arctic.  In 
the  summer,  the  mixed  layer  is  much  deeper  than 
Level  1,  probably  reaching  down  to  Level  4  or  5. 


10.0  LEVITUS  CLIMATOLOGICAL  DATA 

Three-dimensional  salinity  and  temperature  fields 
for  initialization  of  the  PIPS2.0  model  were  prepared 
from  the  Levitus  climatological  data  set.  Seasonal 
and  annual  data  were  used  to  make  monthly 
temperature  and  salinity  estimates,  which  were 
interpolated  into  the  model  domain  according  to  the 
15  levels  of  bottom  topography.  Figures  10.1  and 
10.2  show  the  top  mixed-layer  temperature  through¬ 
out  the  year  and  the  salinity  of  the  mixed  layer, 
respectively,  for  the  PIPS2.0  domain. 

The  interpolation  procedure  is  as  follows.  The 
monthly  data  are  extrapolated  from  the  sea  grid  to 
nearby  land  gridpoints.  The  extrapolated  data  are 
then  used  to  compute  three-dimensional  cubic  spline 
functions  in  longitude,  latitude,  and  depth.  The  cubic 
splines  are  then  used  to  interpolate  the  data  set  onto 
the  PIPS2.0  grid.  The  horizontal  spatial  resolution 
of  the  Levitus  data  is  1°  x  1°  in  the  oceans.  There¬ 
fore,  when  those  data  are  extrapolated  on  to  the 
PIPS2.0  grid  (0.28575°  x  0.28575°)  near  coasts, 
narrow  straits,  channels,  or  semi-enclosed  seas,  the 
accuracy  becomes  questionable. 


A  review  of  the  data  and  general  trends  from  the  Levitus  data  set  may  help  in  understanding 
the  study  results  to  be  presented  in  Sec.  13.0.  Water  temperatures  in  the  Arctic  generally  decrease 
with  increasing  latitude.  The  water  temperature  is  usually  2°-  6°  warmer  in  summer  than  in  winter 
throughout  the  domain,  except  where  permanent  sea  ice  is  present  (the  Central  Arctic).  The  tem¬ 
perature  there  remains  at  the  freezing  point  for  most  of  the  year.  The  salinity  is  3-6  ppt  lower  in 
the  western  Arctic  than  in  the  eastern  Arctic.  The  western  Arctic  obtains  fresh  water  from  the 
Siberian  and  Alaskan  rivers,  while  the  eastern  Arctic  is  mixed  with  salty  North  Atlantic  water. 
The  salinity  in  the  North  Pacific  is  lower  than  in  the  North  Atlantic  Ocean,  the  Norwegian  Sea, 
and  the  East  Greenland  Seas. 


11.0  RIVER  RUNOFF 

The  ice-ocean  coupled  model  includes  runoff  from  eight  major  rivers  (Table  11.1),  as  suggested 
by  Aagaard  and  Carmack  [1989].  The  Severnaya  Dvina  River  is  not  used  in  the  modeling  because 
the  White  Sea  to  which  it  flows  is  closed  from  other  Arctic  seas.  Based  on  Shiklomanov  and 
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Fig.  10.1  — Monthly  Levitus  climatological  data  of  temperature  (‘’C)  at  15  m 
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Fig.  10.1 — cont. 
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Fig.  10.2  —  Monthly  Levitus  climatological  data  of  salinity  (ppt)  at  15  m 
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Fig.  10.2  —  com. 
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River 

Runoff 

(km^yr-*) 

MacKenzie 

340 

Kolyma 

102 

Indigirka 

57 

Lena 

520 

Kotuy 

105 

Yenisei 

603 

Ob 

530 

Pechora 

130 

Table  11.1  — Major  River  Runoff  in  the  Arctic  [Aagaard  and 

Carmack  1989] 


Skakalsky  [1993],  the  annual  runoff  from  each  river  was  interpolated  into  monthly  data  to  be 
consistent  with  the  other  data  fields  (the  NOGAPS  winds  and  the  Levitus  climatology)  that  were 
used  for  studies.  Flood  or  drought  seasons  would  change  the  monthly  river  flow,  which  in  turn 
would  change  both  ice  and  ocean  conditions  near  river  mouths. 


12.0  COUPLING  THE  ICE  AND  OCEAN  MODELS 

The  discussion  to  this  point  has  left  the  interface  between  the  Cox  ocean  model  and  the  ice 
growth  model  undefined.  In  the  earlier  versions  of  the  PIPS  model  [Preller  1985],  the  ice  model 
essentially  ran  independently  of  the  ocean  model.  The  oceanic  forcing  necessary  for  the  former 
model  was  provided  by  Hibler  and  Bryan  [1987].  Changes  in  salinity,  the  daily  exchanges  of  heat 
between  ice  and  ocean,  and  the  drag  that  the  ice  exerts  on  the  currents  directly  under  it  were  totally 
unaccounted  for  in  earlier  work.  The  ice-ocean  coupled  model  described  here  has  sought  to  improve 
that  situation  by  using  an  algorithm  that  simultaneously  integrates  both  the  ocean  model  and  the  ice 
model,  while  exchanging  information  along  the  way.  That  algorithm  is  described  in  this  section. 

The  intent  of  the  algorithm  is  to  preserve  the  autonomy  of  both  the  ice  model  and  the  Cox 
ocean  model,  and  couple  the  models  by  exchanging  necessary  information.  The  top  mixed-layer 
temperatures,  the  variable  freezing  temperatures,  the  oceanic  heat  fluxes,  and  the  ocean  currents  are 
provided  from  the  ocean  model  to  the  ice  model.  In  return,  a  variable  drag  coefficient,  salinity 
transfers,  and  daily  ice  growth  rates  are  provided  to  the  ocean  model.  These  interchanges  are  made 
on  each  day  of  model  time.  The  details  of  these  transfers  follow. 

The  two  important  equations  in  the  ocean  model  into  which  the  daily  information  from  the  ice 
model  is  updated  are  the  equations  for  water  temperature  and  salinity.  Following  Cheng  and  Preller 
[1992;  1994].  Hibler  and  Bryan  [1987],  Cox  [1984],  and  Sarmiento  and  Bryan  [1982],  the  heat 
equation  for  computing  water  temperature  is  written  as  follows: 


dT 


dz^ 


-RtiT-To), 


^mix 


(4) 


where  T  is  the  water  temperature, 
u  is  the  ocean  current, 

Kfi  is  the  coefficient  of  the  horizontal  eddy  diffusion. 
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is  the  horizontal  Laplacian  operator, 

Kr  is  the  coefficient  of  the  vertical  eddy  diffusion, 

fA  is  the  ice  growth/melting  rate  in  open  water  due  to  atmospheric  forcing  only, 

A  is  the  ice  concentration, 

5(z)  is  the  delta  function  (i.e.,  one  in  the  mixed-layer  and  zero  otherwise), 

Rq  is  a  ratio  of  the  latent  heat  of  fusion  of  sea  ice  to  the  heat  capacity  of  seawater, 

Q{T-  Tj)  is  1  when  the  mixed-layer  temperature  T  is  greater  than  the  freezing  point  Tf,  and 
zero  otherwise, 

'  is  the  robust  constraint  for  water  temperature  and  salinity,  and  is  250  days  for  all  levels 
in  the  ice-ocean  coupled  model. 

To  is  the  Levitus  monthly  climatological  temperature,  and 

Zmix  is  the  mixed-layer  thickness. 

The  robust  constraint,  R^  * ,  was  used  in  the  integration  of  the  equation,  not  only  to  keep  ocean 
temperature  and  salinity  from  dissipating  due  to  eddy  diffusion  and  lack  of  precipitation,  but  also 
to  allow  atmospheric  heating/cooling  to  penetrate  into  the  mixed  layer  in  open  water  areas.  This 
constraint  differs  from  those  used  in  the  Cox  ocean  model  (Sec.  8.0),  which  were  imposed  to 
accelerate  the  model  spin-up  from  a  motionless  ocean,  as  well  as  to  relax  the  temperature  and 
salinity  towards  the  Levitus  climatological  data.  The  sign  convention  used  here  is  that  positive  ice 
growth  rate  due  to  atmospheric  forcing  fA  means  that  the  water  is  cooling  off;  negative  fA  implies 
the  water  is  warming. 


The  salinity  equation  in  the  ocean  model  is: 


as 


a2s  0.035S.6(z) 


-R,(S-So), 


(5) 


where  S  is  the  salinity, 

Sf  is  the  total  ice  growth  rate  in  open  water,  and 
So  is  the  Levitus  monthly  climatology  salinity. 

A  robust  constraint  is  also  used  here  for  the  same  purpose  as  with  the  temperature  equation. 
Salinity  in  the  mixed  layer  could  vary  seasonally  according  to  ice  growth  or  melt  because  sea  ice 
rejects  salt  as  it  forms.  An  average  salinity  value  of  35  ppt  is  used  for  fresh  water  when  melting 
occurs  and  for  brine  when  ice  is  growing. 

For  numerical  stability  of  the  ocean  model  when  it  is  coupled  with  the  ice  model,  a  0.5-hr 
timestep  is  used  for  integrating  the  temperature  and  salinity  equations,  and  a  0.05-hr  timestep  is 
used  for  the  stream  function  and  ocean  current.  For  the  ice  model,  a  2-hr  timestep  was  chosen. 
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The  freezing  temperature  (°C)  depends  on  the  mixed-layer  salinity,  as  follows. 

Tf=  -54.4  S,  where  S  is  the  salinity. 

The  mixed-layer  oceanic  heat  fluxes  are  defined  as  the  heat  advected  and  diffused  into  each 

2  d^T 

grid  cell  and  have  the  following  forms:  -  V'(ur),  and/sT,— ,  each  multiplied  by  the 

water  mass  and  its  heat  capacity.  Vertical  mixing  and  vertical  heat  convection  in  the  Cox  ocean 
model  takes  place  whenever  density  differences  occur  that  must  be  smoothed  out.  In  this  coupled 
ice-ocean  model,  mixed-layer  currents  that  are  driven  by  geostrophic  winds  are  used  instead  of  the 
second-layer  ocean  currents  [Hibler  and  Bryan  1987;  Semtner  1987]. 

In  this  study,  a  variable  drag  coefficient  between  sea  ice  and  water  was  applied  according  to 
the  boundary  layer  theory  of  McPhee  [1990]: 


Cd,ice-water  — 


where  Oddce-watec  ts  the  coefficient, 
k  is  the  von  Karmon  constant,  0.41, 
zo  is  the  ice  roughness  and  is  set  to  be  0.01  m,  and 
h  is  the  ice  thickness  and  is  set  to  be  greater  than  zq  times. 

When  h  is  equal  to  2.5  m,  the  drag  coefficient  is  equal  to  the  constant  value  (0.0055)  used  by 
the  initial  versions  of  PIPS.  Note  that  the  drag  coefficient  decreases  as  the  ice  thickness  increases; 
i.e.,  there  is  very  little  drag  stress  on  thick  ice  from  water.  However,  the  drag  stress  becomes  the 
dominant  factor  over  wind  stress  when  h  is  small  and  is  very  close  to  zo- 

Hibler  and  Bryan  [1987]  used  the  ice  internal  stress  and  the  surface  wind  stress  in  the  momentum 
balance  equation  of  the  ocean  model.  In  coupling  the  PIPS2.0  ice  and  ocean  models,  the  sea  ice  is 
treated  as  a  boundary  layer  that  blocks  direct  heat  and  momentum  exchanges  between  the  atmosphere 
and  the  ocean.  Surface  wind  stress  passes  momentum  to  the  sea  ice;  some  momentum  moves  the 
sea  ice,  and  some  transfers  into  the  internal  ice.  Sea-ice  motion  that  drags  the  top  mixed  layer  at 
the  surface  is  modeled  in  PIPS2.0. 

In  the  coupled  model,  the  mixed-layer  temperature  depends  on  the  open-water  condition.  If 
the  open  water  starts  to  grow  sea  ice,  then  the  mixed-layer  temperature  is  set  at  the  freezing 
point.  The  open  water  does  not  grow  sea  ice  if  the  atmospheric  and  oceanic  forcings  are  not  cold 
enough.  They  simply  reduce  the  mixed-layer  temperature.  In  the  melting  seasons,  the  forcings 
warm  the  mixed  layer  or  melt  the  existing  sea  ice. 

Since  the  Hibler  and  PIPS  ice  models  do  not  couple  with  ocean  models,  the  mixed-layer 
temperature  depends  on  accumulated  heat  during  melting  seasons,  called  YNEG  in  the  model 
Fortran  code.  In  modeling,  the  heat  is  stored  and  accumulated.  It  has  to  be  balanced  completely  in 
fall  and  winter  before  growing  any  sea  ice.  In  the  coupled  model,  the  oceanic  heat  fluxes  were 
computed  based  on  the  heat  advection  and  diffusion  (Eq.  4).  In  the  ocean  model,  the  mixed-layer 
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fAHz)RoQ(T-Tf) 

temperature  is  heated  or  cooled  according  to  the  atmospheric  forcing,  i.e., - — - —  .  Therefore, 

it  is  reasonable  to  eliminate  the  YNEG  term  in  the  coupled  model. 

Figures  12.1-12.6  show  the  ice-ocean  coupled  model  in  action.  The  monthly  Levitus  data  and 
the  monthly  NOGAPS  wind  forcing  data  were  used  to  drive  the  ice-ocean  coupled  model  for  5  yr. 
The  model  continued  to  run  for  another  2  yr  with  the  monthly  river  runoff  included.  Figures  12.1 
and  12.2  show  the  model  distributions  of  ice  thickness  and  concentration  for  each  month  of  a  year. 
The  model  ice-drift  velocities  for  the  same  year  are  plotted  in  Fig.  12.3.  Figures  12.4  and  12.5 
are  the  model  distributions  of  ocean  temperature  and  salinity,  and  Fig.  12.6  is  the  ocean  current  at 
a  15-m  depth. 

For  comparison,  the  dark  thick  line  in  the  March  ice  concentration  of  Fig.  12.2  indicates 
ice-edge  locations  near  mid-month  that  were  derived  as  the  weekly  ice  concentration  analysis  by 
the  Navy-NOAA  (National  Oceanic  and  Atmospheric  Administration)  Joint  Ice  Center  (JIC).  The 
JIC  weekly  analysis  is  derived  from  satellite  passive  microwave  measurements,  visible  and  infrared 
imagery,  aerial  reconnaissance  and  conventional  field  observations.  The  visual  agreement  between 
the  data  and  the  model  results  is  good. 


13.0  RESULTS  OF  AN  ICE-OCEAN  COUPLED  MODEL 

Section  13.1  shows  the  ocean  model  results  of  the  PIPS2.0  ice-coupled  model  and  discusses 
how  the  parameters  of  viscosity  and  diffusivity,  the  lateral  wall  boundaries,  and  the  forcing  fields 
influence  ocean  circulation.  After  that,  for  further  regional  discussion,  the  PIPS2.0  model  domain 
was  divided  into  seven  regions  in  the  area  of  interest  as  follows;  the  Central  Arctic  (Sec.  13.2),  the 
Barents  Sea  (Sec.  13.3),  the  Norwegian/East  Greenland  Seas  (Sec.  13.4),  Baffin  Bay/Labrador  Sea 
(Sec.  13.5),  Hudson  Bay  (Sec.  13.6),  the  Bering  Sea  (Sec.  13.7),  and  the  Sea  of  Okhotsk  (Sec.  13.8). 
Since  each  of  these  regions  has  the  potential  of  becoming  a  separate  topic  of  interest  in  future 
international  affairs,  each  is  discussed  separately,  and  the  best  available  model  description  is  given. 


13.1  The  Ocean  Model  Results 

In  the  Arctic  and  sub-Arctic  regions,  all  the  major  circulations  appear  in  the  output  of  the 
ice-ocean  coupled  model:  the  Alaska  Gyre,  the  Oyashio  Current,  the  Beaufort  Gyre,  the  Transpolar 
Drift,  the  East  Greenland  Current,  the  Norwegian  Current  and  its  branches  to  the  Barents  Sea  and 
south  and  west  of  Spitsbergen,  the  North  Atlantic  Drift  and  its  return  current,  and  the  Labrador 
Current.  In  general,  the  stream  function  follows  the  f/H  contours  throughout  the  domain,  as  Semtner 
[1987]  suggested,  where  /  is  the  Coriolis  parameter  and  H  is  the  ocean-bottom  topography  (see 
Figs.  13.1.1  and  9.1).  Figure  13.1.1  shows  the  stream  functions  in  March,  June,  September,  and 
December.  The  features  in  the  stream  functions  remain  almost  constant  through  the  year,  as  expected, 
because  the  values  represent  the  major  ocean  circulation.  The  seasonal  variations  represent  the 
minor  differences  between  the  features. 

The  original  values  of  the  horizontal  coefficients  of  eddy  viscosity  and  diffusivity,  which  were 
used  in  the  Cox  ocean  model,  were  appropriate  when  the  model  was  run  at  2°-3°  spatial  resolution. 
These  coefficients  were  quite  large,  however,  for  the  fine  resolution  (0.28575°)  of  the  coupled 
model.  Sensitivity  studies  to  determine  better  coefficients  for  both  ice  and  ocean  models  will  be 
conducted  and  results  presented  in  the  future. 
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Fig.  12.1  — Model  ice  thickness  (m)  distribution,  Jan-Dec  1986 
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Fig.  12,1 — cont. 
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Fig.  12.2  —  Model  ice  concentration  distribution,  Jan-Dec  1986 
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Fig.  12.2  —  cont. 
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Fig.  12.2  —  cont. 
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Fig.  12.3  —  Model  ice  drift  velocity  (cm/s),  Jan-Dec  1986 
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Fig.  12.3  —  cont. 
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Fig.  12.3  —  cont. 
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Fig.  12.4  —  Model  ocean  temperature  (°C)  at  15  m,  Jan-Dec  1986 
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Fig.  12.4  —  cont. 
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Fig.  12.5  —  Model  ocean  salinity  (ppt)  at  15  m,  Jan-Dec  1986 
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Fig.  12.6  —  Model  ocean  currents  (cm/s)  at  15  m,  Jan-Dec  1986 


The  Development  of  an  Ice-Ocean  Coupled  Model 


0.200  Ef02 
MAXIMUM  VECTOR 


198603 


198608 


Fig.  12.6  —  cont. 
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Fig.  12.6  —  cont. 
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Fig.  13.1.1  — Model  stream  functions  (cm“/s)  (a)  Mar,  (b)  Jun,  (c)  Sep,  and  (d)  Dec  1986 
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Because  of  the  way  this  model  is  implemented,  inaccuracies  are  introduced  into  some  of  the 
calculations  at  lower  latitudes.  For  example,  inaccurate  ocean  currents  on  the  lower  right,  upper  left 
and  lower  left  corners  of  Fig.  12.6  (i.e.,  near  the  Gulf  Stream,  the  Kuroshio  Current,  and  the 
California  Current,  respectively)  are  caused  by  (1)  an  assumption  in  the  model  that  operates  mainly 
in  those  lower  latitudes,  (2)  the  lateral  solid  wall  boundary,  and  (3)  the  mass  conservation  law.  The 
solid  wall  assumption  blocks  water  exchanges  between  high-  and  low-latitude  water  masses.  As  a 
result,  the  model  California  Current  moves  northward,  opposite  to  its  observed  motion.  However, 
since  the  boundary  at  which  the  solid  wall  assumption  is  needed  has  been  placed  far  away  from  sea 
ice-covered  regions,  the  influence  of  the  assumption  on  ice  growth  and  drift  should  be  negligible. 

Another  check  on  the  current  predictions  of  the  PIPS2.0  model  has  recently  become  available, 
but  again,  that  check  is  at  low  latitudes.  According  to  Ebbesmeyer  and  Ingraham  [1994],  ocean 
current  pathways  tracked  by  following  the  trajectories  of  toys  and  shoes  spilled  in  a  maritime 
accident  start  from  the  middle  of  the  North  Pacific,  move  eastward  to  the  coast  of  North  America, 
and  finally  flow  northward  to  Alaska.  Another  route  starts  at  a  low-latitude  location  (close  to 
135°  W,  13°  N),  moves  eastward  along  that  latitude,  then  turns  northeast  toward  Baja  California  (at 
(115°  W,  15°  N),  and  finally  travels  northward  along  the  California  coast  to  Alaska.  Although  this 
route  is  consistent  with  the  PIPS2.0  model  current  pathways,  that  could  be  a  special  case  driven  by 
unusual  winds  or  coastal  Kelvin  waves. 

The  current  ocean  model  in  PIPS2.0  cannot  be  used  for  resolving  eddies.  Because  of  the  high 
horizontal  eddy  coefficients  just  mentioned  and  the  model’s  horizontal  grid  resolution,  seasonal 
variations  of  regional  meanders  and  eddies,  as  well  as  those  caused  by  ice  melt  and  growth,  cannot 
be  resolved  with  the  PIPS2.0  model.  For  example,  the  Chukchi  Sea  front  between  the  Bering  Sea 
warm  inflow  and  the  Arctic  cold  water  cannot  be  delineated  clearly  by  looking  at  Figs.  12.4,  12.5, 
and  12.6,  even  though  the  model  temperature  and  salinity  vary  seasonally.  The  noneddy-resolving 
Cox  ocean  model  is  used  to  compute  large-scale  ocean  circulations  and  the  water  transport  between 
high  and  low  latitudes. 

.  The  results  of  all  these  studies  are  very  sensitive  to  variations  in  the  wind  fields.  The  model 
ocean  currents,  when  driven  by  the  monthly  winds,  change  slowly  during  the  year.  When  the 
NOGAPS  daily  (3-  or  6-hourly)  winds  are  used,  however,  the  currents  show  strong  daily  variations, 
especially  in  the  top  mixed  layer.  For  example.  Fig.  13.1.2  shows  such  ocean  currents  and  salinities 
for  21  and  31  October  1992,  when  the  coupled  model  was  driven  by  6-hourly  wind  forcing  fields. 
Wind-driven  currents  in  open  water  are  more  energetic  and  variable  than  the  monthly  wind  case, 
even  when  the  salinity  distributions  are  the  same.  More  sensitivity  studies  of  the  wind  stresses  and 
the  air  drag  coefficient  will  be  conducted  in  the  future. 

Although  the  mixed-layer  temperature  starts  with  the  Levitus  climatology,  a  loose  constraint  is 
applied  during  the  run  integration.  This  constraint  allows  the  atmospheric  forcing  to  significantly 
change  the  temperature  with  season.  For  example,  the  Norwegian  Sea  reaches  8°C  in  March  and 
1°C  in  September.  Similarly,  the  mixed- layer  salinity  also  shows  seasonal  variations.  In  the  East 
Siberian  Sea  and  the  Beaufort  Sea,  the  model  salinity  is  about  1-2  ppt  higher  than  the  Levitus  data. 
Seasonal  variations  of  ocean  currents  and  temperature  and  salinity  at  depth  are  much  smaller  than 
those  in  the  mixed  layer. 

Meunch  [1989]  classified  the  Arctic  Ocean  and  the  Norwegian  and  East  Greenland  Seas  into 
four  different  categories  of  water:  the  Arctic  surface  layer  (ASL),  the  Atlantic  layer,  the  deep 
waters,  and  the  Arctic  shelf  waters.  The  three  main  classes  of  ASL  can  be  recognized  as  follows: 
the  Atlantic  Water  (AW),  which  is  the  Norwegian-Atlantic  Current  with  temperatures  above  3°C 
and  salinity  above  34.9  ppt;  the  Polar  Water  (PW),  which  has  been  diluted  by  admixtures  with  fresh 
water,  so  its  temperature  and  salinity  are  below  0°C  and  34.4  ppt,  respectively;  and  the  Arctic 
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Fig.  13.1.2  —  Model  ocean  temperature  (°C)  and  current  (cm/s)  on  21  and  31  Oct  1992 
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Surface  Water  (ASW),  which  has  temperatures  between  0°-  3°C  and  salinity  between  34.4-34.9  ppt. 
AW  is  located  basically  along  the  Norwegian  Current  and  PW  in  the  central  Arctic  Ocean  and  the 
East  Greenland  Current.  ASW  lies  between  AW  and  PW  and  in  the  gyres  of  the  Greenland  and 
Iceland  Seas. 

The  model  distributions  of  temperature  and  salinity  (Figs.  12.4  and  12.5)  agree  well  with  the 
three  ASL  water  types.  Note  that  the  model  temperature  and  salinity  vary  seasonally.  For  example, 
the  summer  temperature  of  the  East  Greenland  Sea  (Fig.  12.4)  is  1°-2°C  higher  than  the  PW 
described  above.  The  model  Atlantic  layer,  deep  waters,  and  Arctic  shelf  waters  are  located  under 
the  mixed  layer  and  will  be  discussed  with  other  sensitivity  studies  in  the  future. 


13.2  Central  Arctic 

The  Transpolar  Drift  and  the  Beaufort  Gyre  are  the  dominant  ocean  currents  in  this  region.  The 
vertical  and  lateral  heat  fluxes  are  the  important  factors  as  to  whether  or  not  sea  ice  exists.  Because 
of  the  permanent  existence  of  sea  ice,  the  mixed-layer  temperature  of  the  currents  remains  at  the 
freezing  point  most  of  the  year.  Consequently,  lateral  oceanic  heat  fluxes  are  small.  However, 
coastal  seas  along  Russia  and  Alaska  are  clear  of  sea  ice  in  summer  and  have  open  water. 

The  distribution  of  the  model  ice  thickness  shows  packed  thick  sea  ice  north  of  Canada  and 
Greenland  (between  3.5  -  8.5  m  for  March  1986).  The  sea  ice  becomes  thinner  toward  the  North 
Pole  (between  3.0  -  4.5  m).  The  thickness  continues  to  decrease  toward  the  Laptev  Sea  and  the 
Kara  Sea  north  of  Russia.  The  distribution  remains  approximately  the  same  throughout  the  year,  but 
the  absolute  value  of  the  thickness  varies  seasonally  (e.g.,  1-2  m  between  summer  and  winter  near 
the  North  Pole). 

The  ice  drift  is  essentially  dominated  by  the  geostrophic  winds  that  blow  from  Russia  to 
Canada  most  of  the  year.  Such  ice  advection  affects  the  ice  distribution;  thick  near  Greenland  and 
thin  near  the  Laptev  Sea.  The  accumulated  thick  ice  in  the  East  Siberian  Sea  is  likely  to  have  been 
caused  by  the  use  of  the  monthly  mean  winds  advecting  the  sea  ice  in  that  direction  for  several 
months.  The  monthly  geostrophic  winds  of  summer  move  from  the  central  Arctic  to  the  Laptev  Sea 
and  force  ice  to  drift  into  that  sea.  As  a  result,  the  Laptev  Sea  is  covered  by  sea  ice  during  summer 
(Figs.  12.1,  12.2,  and  12.3). 

In  similar  work,  Walsh  et  al.  [1985]  used  the  Hibler  ice  model  to  extend  the  model  domain  to 
the  North  Pacific  Ocean  and  the  Labrador  Sea.  The  horizontal  grid  was  222  km  and  the  model  sea 
ice  was  estimated  in  Cartesian  coordinates.  No  ocean  model  was  coupled.  Geostrophic  winds  and 
temperature-derived  thermodynamic  fluxes  from  1951  to  1980  drove  the  model.  The  30-yr  mean  for 
February  (Case  N30)  shows  an  ice  thickness  distribution  similar  to  the  PIPS2.0  result:  thin  in  the 
Laptev  Sea  (approximately  1.5  m),  thick  near  the  North  Pole  (approximately  2.5  m),  and  still 
thicker  toward  Greenland  and  the  Canadian  Archipelago  (approximately  3.5-4.0  m).  The  August 
30-yr  mean  shows  a  similar  distribution,  but  thinner:  0-0.5  m  for  the  Laptev  Sea,  2.5-3.0  m  north 
of  Greenland  and  the  Canadian  Archipelago.  In  spite  of  different  atmospheric  forcing  fields,  and 
without  the  ocean  model,  Walsh  et  al.’s  February  and  August  means  are  consistent  with  the  PIPS2.0 
ice-ocean  coupled  model  results,  although  the  former  is  0.5-4.0  m  thinner  than  the  latter. 


13.3  Barents  Sea 

The  ice  thickness  reaches  a  maximum  of  2  m  in  the  northwestern  Barents  Sea  in  March,  and 
the  ice  drift  is  dominated  by  the  cyclonic  geostrophic  winds.  No  sea  ice  grows  in  the  south 
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and  southwestern  Barents  Sea  because  the  water,  which  is  part  of  the  Norwegian  Current,  is  warm  and 
salty  throughout  the  year.  In  March,  the  open-water  temperature  ranges  from  0°-4°C,  even  though 
the  corresponding  air  temperature  is  between  —10°  and  0°C.  In  summer,  sea  ice  usually  retreats 
from  the  Barents  Sea  to  a  region  near  the  Franz  Josef  Islands,  the  boundary  of  the  sea  and  the 
central  Arctic.  The  mixed  layer  is  dominated  by  Atlantic  water  and  has  a  salinity  between  34  and 
35  ppt,  except  for  the  northern  sea  near  Spitsbergen,  the  Franz  Josef  Islands,  and  Novaya  Zemlya. 
That  region  could  have  a  salinity  as  low  as  33  ppt  in  summer.  The  sea  ice  distribution  affected  by 
geostrophic  winds  could  vary  from  year  to  year. 

In  the  Barents  Sea,  the  ice  edge  of  the  PIPS2.0  model  agrees  in  general  with  the  JIC  weekly 
analysis  throughout  the  year.  This  consistency  implies  that  the  ice-ocean  coupled  model  gives  the 
proper  currents  and  heat  fluxes,  especially  near  the  ice  edge.  The  JIC  analysis  reveals  two 
inconsistencies.  For  August  and  September,  the  model  results  show  extra  sea  ice  east  of  Spitsbergen. 
For  October  and  November,  the  model  results  indicate  fast  ice  growth  (generally  less  than  20%  in 
concentration)  in  the  central  and  northern  portions  of  the  sea.  These  inconsistencies  might  imply 
that  the  monthly  atmospheric  forcing  does  not  reflect  such  sea  ice  change  in  fast  melting  and 
growing  seasons.  When  the  3-hourly  atmospheric  forcing  is  used  in  operation,  the  model  sea  ice 
would  follow  the  temporal  variations  and  the  inconsistencies  should  be  improved. 


13.4  East  Greenland  and  Norwegian  Seas 

The  East  Greenland  Current  brings  cold,  fresh  Arctic  surface  water  into  the  seas  through  the 
Fram  Strait.  This  current  travels  farther  south,  passes  the  Denmark  Strait,  and  enters  the  North 
Atlantic  Ocean.  The  water  is  less  saline  in  summer  due  to  ice  melt  in  the  sea  and  the  Arctic  Ocean 
(approximately  33  ppt),  but  is  saltier  in  winter  due  to  ice  growth  (approximately  34  ppt).  The 
Norwegian  Current,  which  is  part  of  the  salty,  warm  Atlantic  water,  flows  from  the  Faroe-Shetland 
Channel  north  to  Spitsbergen  and  remains  ice-free  during  the  year. 

The  strong  ocean  current  and  wind  off  the  eastern  coast  of  Greenland  move  sea  ice  from 
the  Fram  Strait  through  the  Denmark  Strait  to  the  southern  tip  of  Greenland  (i.e.,  from  about 
82°  N-60°  N).  In  winter,  sea  ice  covers  the  entire  region;  in  summer,  it  retreats  as  far  south  as 
the  Denmark  Strait.  In  March,  for  example,  the  water  temperature  north  of  the  Denmark  Strait 
can  reach  4°C.  However,  due  to  the  strong  advection  in  this  region,  sea  ice  can  still  survive. 

In  August,  the  southwesterly  offshore  winds  near  the  Fram  Strait  advect  sea  ice  away  from  the 
coast.  As  a  result,  the  distribution  of  model  sea  ice  indicates  open  water  there  (or  implies  possible 
fast  ice).  The  JIC  weekly  analysis  shows  fast  sea  ice,  but  no  open  water  along  the  coast.  The  model 
open  water  is  probably  the  result  of  using  the  offshore  winds  to  drive  the  coupled  model  for  the 
entire  month  of  August,  while  the  JIC  observations  are  taken  on  only  a  few  days.  That  sea  ice 
distribution  from  the  model  might  be  different  if  the  model  was  driven  by  daily  (or  3-hourly)  winds, 
instead  of  the  monthly  winds.  If  daily  winds  were  used,  sea  ice  might  not  be  constantly  expelled 
from  the  coast.  Open  water  in  the  model  does  not  exist  when  using  the  1992  daily  winds  for  August 
(P.  Posey,  pers.  comm.). 


13.5  Baffin  Bay  and  Labrador  Sea 

After  picking  up  some  Atlantic  water  southwest  of  Iceland,  the  East  Greenland  Current  flows 
around  the  southern  tip  of  Greenland  and  into  the  Labrador  Sea.  It  is  renamed  as  the  West  Greenland 
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Current  (WGC)  and  continues  northward  along  the  west  coast  of  Greenland.  Most  of  the  WGC 
water  drifts  to  the  west  and  mixes  with  the  Baffin  Land  Current  and  the  Labrador  Current.  The 
WGC  contains  warm  Atlantic  water  and  remains  ice-free  along  its  path  throughout  the  year.  In 
summer,  the  water  temperatures  of  Baffin  Bay  and  the  Labrador  Sea  seem  to  decrease  with  the 
distance  from  the  North  Atlantic  Ocean.  The  water  salinity  ranges  between  33-35  ppt  in  winter  and 
between  32-34  ppt  in  summer. 

Baffin  Bay  and  the  Labrador  Sea  have  very  little  sea  ice  in  summer  and  grow  new  sea  ice  every 
winter.  Advected  sea  ice  from  the  Arctic  Ocean  is  limited  because  of  the  narrow  straits,  channels, 
and  sounds.  In  March,  Baffin  Bay  is  covered  by  sea  ice,  except  for  areas  near  the  west  Greenland 
coast  where  offshore  winds  constantly  advect  sea  ice  away.  As  with  the  model  predictions  for  the 
East  Greenland  and  Norwegian  Seas,  when  the  1992  daily  winds  are  used  in  studies,  sea  ice  entirely 
covers  Baffin  Bay  (P.  Posey,  pers.  comm.).  The  model  ice  edge  agrees  with  the  JIC  analysis  except 
in  two  regions:  Baffin  Bay  west  of  Greenland  and  the  northern  Labrador  Sea  south  of  the  Ellesmere 
Island.  The  absence  of  sea  ice  in  the  model  results  is  likely  to  have  been  caused  by  the  constant 
offshore  monthly  winds,  as  discussed  above.  According  to  the  JIC  data  in  the  southern  Labrador 
Sea,  sea  ice  should  exist  along  the  Canada  coast  and  extend  to  Newfoundland.  In  the  coupled 
model,  March  cyclonic  winds  move  young  sea  ice  from  north  to  south  and  then  from  the  southern 
Labrador  Sea  into  the  warm  Atlantic  water. 

Comparing  the  March  ice  thickness  with  the  March  ice  concentration  reveals  that  the  latter  ice 
edge  extends  farther  toward  the  North  Atlantic  Ocean  than  the  former.  The  difference  in  the  two 
ice-edge  locations  implies  thin  sea  ice  in  the  southern  Labrador  Sea.  Thin  sea  ice  might  be  caused 
by  returned  warm  Atlantic  water,  which  plays  an  important  role  in  melting  sea  ice. 


13.6  Hudson  Bay 

In  summer,  the  model  water  temperature  in  Hudson  Bay  varies  from  0°-10°C  and  its  salinity 
ranges  from  28-32  ppt.  Seasonal  river  runoff,  which  contributes  fresh  water  and  causes  low  salinity  in 
the  east  and  south  of  the  bay,  was  not  included  in  this  modeling.  River  runoff  could  affect  the 
horizontal  stratification  of  water  in  the  bay  and  produce  an  estuarine-type  circulation  [Pickard  and 
Emery  1990].  In  winter,  water  salinity  increases  to  30-33  ppt,  mainly  due  to  ice  growth  and  brine 
rejection  into  the  bay. 

Hudson  Bay  is  totally  covered  by  sea  ice  during  winter  and  is  completely  free  of  ice  in  August 
and  September.  The  bay  is  partially  occupied  by  sea  ice  in  spring  and  fall.  Differences  in  the 
estimates  of  open  water  between  the  ice  thickness  and  ice  concentration  distributions  represent  thin 
sea  ice  coverage,  especially  in  November  and  December  when  it  grows  rapidly.  Hudson  Bay  receives 
little  influence  from  warm  North  Atlantic  water,  except  for  the  northern  part  of  the  bay,  where  some 
inflow  from  the  Labrador  Sea  can  occur  [Wang  et  al.  1994a;  1994b].  Most  of  the  ice  growth  and 
melt  in  Hudson  Bay  are  dominated  by  atmospheric  forcing. 


13.7  Bering  Sea 

The  Bering  Sea  inflow  comes  from  the  Aleutian  Islands  close  to  the  Alaska  Peninsula.  The 
water  turns  westward  and  flows  out  of  the  sea  through  the  western  Aleutians.  This  outflow  joins 
the  Oyashio  Current  east  of  the  Kamchatka  Peninsula.  The  mixed-layer  temperature,  influenced  by 
North  Pacific  water,  is  sometimes  higher  than  6°C  in  the  summer.  The  water  exchange  with  the 
Arctic  through  the  Bering  Strait  does  not  seem  to  have  significant  impact  for  the  sea  when  com¬ 
pared  with  that  of  the  North  Pacific  Ocean,  but  it  could  possibly  initiate  and  then  accelerate  sea  ice 
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growth  near  the  Strait.  Little  sea  ice  is  advected  from  the  Arctic  through  the  Bering  Strait.  Similar 
to  the  situation  with  temperature,  the  mixed-layer  salinity  depends  mostly  on  North  Pacific  inflow 
and  brine  rejection  due  to  sea  ice  growth. 

In  the  Bering  Sea,  sea  ice  starts  to  form  in  late  October  or  early  November  and  starts  to  melt 
in  April.  The  sea  remains  open  throughout  the  summer.  The  model  ice-edge  locations  agree  well 
with  the  JIC  weekly  analysis,  except  for  the  southeast  Bering  Sea  in  winter.  In  March,  for  example, 
the  model  results  indicate  the  existence  of  sea  ice  in  the  northern  Bering  Sea  and  along  the  western 
Siberian  coast.  Sea  ice  that  drifts  westward  is  basically  dominated  by  easterly  geostrophic  winds. 
In  the  southeast  Bering  Sea,  easterly  offshore  winds  blow  from  Alaska  to  the  sea  and  advect  sea 
ice  away  from  the  coast.  In  the  model,  consequently,  the  southeast  Bering  Sea  cannot  retain  sea  ice. 
However,  the  JIC  weekly  analysis  of  ice  concentration  shows  some  sea  ice  in  this  region,  which 
disagrees  with  the  model  results.  When  variable  1992  daily  winds,  instead  of  the  monthly  winds, 
are  used  to  drive  this  model,  most  of  the  southeastern  sea  is  covered  by  sea  ice  in  winter  [P.  Posey, 
pers.  comm.].  Another  reason  the  model  incorrectly  captures  this  area  is  that  fresh  river  runoff 
makes  sea  ice  grow  easier  because  the  lack  of  salinity  causes  the  freezing  temperature  to  be  higher 
than  in  regular  seawater.  River  runoff  was  not  included  in  this  part  of  the  modeling,  although  there 
is  actually  considerable  outflow.  This  situation  will  be  addressed  in  future  work. 


13.8  Sea  of  Okhotsk 

Some  of  the  Oyashio  Current  has  been  observed  to  make  the  turn  around  the  southern  tip  of 
the  Kamchatka  Peninsula  and  to  flow  into  the  Sea  of  Okhotsk.  It  continues  northward  along  the 
west  coast  of  that  peninsula  and  circulates  the  sea  cyclonically.  Finally,  the  water  flows  out  through 
the  southern  Kuril  Islands.  Another  small  inflow  (1-3  Sv)  comes  in  from  the  Soya  Channel  between 
the  Hokkaido  and  Sakhalin  Islands. 


Fig.  13.8.1 — Ocean  current  in  the  Sea  of  Okhotsk  with 
low  eddy  coefficients,  15-m  depth 


The  circulation  in  the  Sea  of  Okhotsk  that 
the  model  produces,  under  normal  values  for 
the  model  parameters,  flows  anticyclonically 
and  is  not  consistent  with  that  discussed 
above.  However,  when  eddy  coefficients  of 
viscosity  and  diffusivity  are  reduced  to  10* 
and  10^  cm^/s,  respectively,  the  cyclonic  cir¬ 
culation  is  reproduced  by  the  model  and  agrees 
with  the  observed  circulation  (Fig.  13.8.1).  Due 
to  plotting  problems.  Fig.  13.8.1  neither  shows 
values  for  the  Kuroshio  and  Oyashio  Currents 
higher  than  5  cm/s,  nor  represents  the  Kuril 
Islands. 

In  the  Sea  of  Okhotsk,  sea  ice  starts  to 
grow  from  late  November  or  early  December, 
and  it  begins  to  melt  in  April.  During  the 
summer,  the  sea  remains  ice-free.  In  winter, 
sea  ice  stays  in  the  northern  and  western  sec¬ 
tions  of  the  sea  (i.e.,  along  the  continental  coast 
and  along  Sakhalin  Island,  respectively).  The 
eastern  Sea  of  Okhotsk  along  the  Kamchatka 
Peninsula  and  along  the  Kuril  Islands  remains 
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ice-free  throughout  the  year.  Differences  between  the  ice  thickness  and  ice  concentration 
distributions  in  the  northern  part  of  the  sea  represent  thin  sea  ice.  The  model  concentration  for 
March  agrees  with  the  JIC  weekly  analysis. 


14.0  REMARKS  AND  CONCLUSIONS 

An  algorithm  was  developed  for  transforming  the  coordinates  of  the  ice  and  ocean  models  that 
make  up  the  PIPS  system  from  Cartesian  coordinates  to  spherical  coordinates.  This  transformation 
was  made  to  avoid  the  numerical  singularity  at  the  North  Pole  and  to  increase  the  numerical 
stability  in  high  latitudes.  The  spherical  coordinates  have  a  new  equator  that  coincides  with  the 
10°E-170°W  great  circle  and  a  new  north  pole  that  is  located  at  the  intersection  of  the  100°  E 
meridian  and  the  true  equator.  Tests  demonstrated  that  the  transformed  PIPS  model  gives  the  same 
results  as  the  Cartesian  coordinates  PIPS  model. 

The  domain  for  the  spherical  coordinates  PIPS  model  was  extended  southward  to  the 
20°  N-30°  N  region,  depending  on  which  part  of  the  grid  is  considered.  PIPS2.0,  as  this  transformed 
and  regridded  model  is  called,  covers  all  sea  ice  in  the  Northern  Hemisphere.  The  horizontal  grid 
resolution  is  0.28575°  x  0.28575°.  The  Cox  ocean  model  has  also  been  regridded  to  the  new  domain. 

The  ice  and  ocean  models  have  been  coupled  for  their  integration  by  exchanging  daily  information. 
Mixed-layer  temperatures,  variable  freezing  temperatures,  oceanic  heat  fluxes,  and  ocean  currents 
go  from  the  ocean  model  to  the  ice  model.  A  variable  drag  coefficient,  salinity  transfers,  and  daily 
ice  growth  rates  are  returned  to  the  ocean  model  from  the  ice  model.  In  coupling,  the  sea  ice  is 
treated  as  a  boundary  layer  that  blocks  direct  heat  and  momentum  exchanges  between  the  atmosphere 
and  ocean. 

The  runs  that  were  made  with  the  ice-ocean  coupled  model  had  much  input  information  in 
common.  The  ET0P05  bathymetry  data  were  interpolated  into  the  PIPS2.0  model  domain  using 
cubic  splines.  After  interpolation  and  integration,  the  bottom  topography  was  divided  into  15  levels, 
as  was  appropriate  for  the  Cox  ocean  model  [1984].  Level  1  (the  top)  thickness  was  fixed  at  30  m. 

The  1986  NOGAPS  atmospheric  forcing  fields  were  used  to  develop  and  test  the  PIPS2.0 
ice-ocean  coupled  model.  These  fields  have  insufficient  resolution,  both  temporally  and  spatially, 
to  support  accurate  studies  in  small  areas.  When  used  operationally,  they  will  be  replaced  by  high- 
resolution,  3-hourly  winds,  and  their  form  will  be  transferred  to  the  ice-ocean  model  as  stresses. 
Both  these  improvements  will  address  limitations  on  the  applicability  of  PIPS2.0  results. 

The  Levitus  seasonal  and  annual  climatological  data  of  water  temperature  and  salinity  [1982] 
were  interpolated  into  the  PIPS2.0  model  domain.  Resolution  of  the  Levitus  data  makes  it  difficult 
to  consider  coasts,  bays,  small  seas,  and  narrow  straits.  New  climatological  data  are  being  prepared 
that  can  address  this  limitation  as  well  and  can  be  used  by  the  PIPS2.0  model  directly. 

In  addition  to  the  limitations  of  applicability  of  PIPS2.0  caused  by  the  input  data  characteristics, 
implementation  algorithm  characteristics  also  might  impact  on  how  successfully  this  model  can  be 
applied  to  given  situations.  The  top  layer  of  the  ice-ocean  coupled  model  that  represents  the  ocean’s 
mixed  layer  has  been  kept  at  a  constant  30-m  thickness  for  this  version.  In  reality,  the  mixed  layer 
is  a  function  of  season  (thin  in  winter  and  thick  in  summer)  and  is  influenced  by  the  influx  of  fresh 
water  from  river  runoff.  In  the  Arctic  Ocean,  the  mixed  layer  is  estimated  to  vary  from  30-50  m 
in  thickness  [Carmack  1990].  The  approximation  that  was  made  could  become  important  in  the 
marginal  ice  zone  because  a  thin  mixed  layer  can  grow  ice  easier  than  a  thick  one.  In  the  North 
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Pacific  and  Atlantic  Oceans,  the  open-water  mixed  layer  varies  seasonally  from  shallow  to  200  m 
deep.  For  both  oceans,  the  top  30-m  level  does  not  represent  the  mixed  layer,  but  it  is  used  in 
modeling  to  interact  with  the  atmospheric  forcing. 

In  the  PIPS2.0  model,  a  variable  ice-water  drag  coefficient  Q,ce-warer  was  used  to  force  certam 
kinds  of  sea  ice  drift  to  follow  the  underlying  ocean  currents  closely.  Of  particular  concern  are  thin 
sea  ice  and  icebergs.  In  reality,  these  ice  forms  drift  along  with  the  currents  like  buoys,  so  it  is 
imperative  that  the  model  represent  that  situation. 

Hibler  and  Bryan  [1987]  used  the  ocean  current  of  Level  2  to  approximate  the  geostrophic 
currents  and  to  compute  the  sea  surface  height.  At  present  in  the  PIPS2.0  model,  the  ocean  currents 
of  Levels  1  and  2  are  similar,  so  the  sea  surface  heights  based  on  both  look  identical,  which  does 
not  affect  sea  ice  computations.  However,  if  the  daily  or  3-hourly  winds  are  used  to  drive  the 
coupled  model,  the  wind-driven  current  of  Level  1  could  significantly  differ  from  that  in  Level  2. 
Further  work  is  needed  on  this  point  to  minimize  the  effect  on  sea  ice  calculations. 
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